# Alex Gazmararian
# agazmararian@gmail.com

source("analysis/visibility/analysis_config.R")
g <- readRDS(here("data", "output", "visibility_analysis.rds"))
g <- subset(g, !is.na(greenbenefit_bin))
source(here("R", "load_functions.R"))

benefits.mean <- round(weighted.mean(g$greenbenefit_bin, g$wt_green_trim) * 100, 0)
print(paste0("Weighted mean of perceived benefits: ", benefits.mean, "%"))
cat(benefits.mean, file = paste0("output/pnas/stats/benefits_mean.txt"))

# Perceived benefits from green projects----

# By whether people see it
p.benefit.vis <- g %>%
  filter(!is.na(greenbenefit)) %>%
  group_by(greenproj_bin, greenbenefit) %>%
  summarize(
    n = n(),
    weighted_count = sum(wt_green_trim)
  ) %>%
  mutate(prop = weighted_count / sum(weighted_count)) %>%
  mutate(greenproj_bin = ifelse(greenproj_bin == 1, "Visible", "Not Visible")) %>%
  ggplot(aes(x = greenbenefit, y = prop, fill = greenproj_bin)) +
  geom_col(position = "dodge") +
  labs(x = NULL, y = "% answered (weighted)", fill = NULL, title = "Perceived project benefits by visibility") +
  geom_text(aes(label = paste0(round(prop*100, 0), "%")), position = position_dodge(width = 0.9), vjust = -0.5, size = plot_text_size) +
  scale_fill_greenproj() +
  theme_classic(base_size = 10) +
  annotate("text", x = Inf, y = Inf,vjust =1, hjust=1.5,  size = 3,color = "grey40", label = paste0("N=",nrow(subset(g, !is.na(greenbenefit))))) +
  scale_y_continuous(expand = c(0, 0), limits = c(0, 1), labels = scales::percent_format(accuracy = 1)) +
  theme(
    legend.position = "inside",
    legend.position.inside = c(.15, .8)
    )
save(p.benefit.vis, file = here("data", "output", "fig_perceived_benefits_vis.rdata"))

# By partisanship
p.benefit.pid <- g %>%
  filter(!is.na(greenbenefit)) %>%
  group_by(pid3, greenbenefit) %>%
  summarize(
    n = n(),
    weighted_count = sum(wt_green_trim)
  ) %>%
  mutate(prop = weighted_count / sum(weighted_count)) %>%
  ggplot(aes(x = greenbenefit, y = prop, fill = pid3)) +
  geom_col(position = "dodge") +
  geom_text(aes(label = paste0(round(prop*100, 0), "%")), position = position_dodge(width = 0.9), vjust = -0.5, size = plot_text_size) +
  labs(x = NULL, y = "% answered (weighted)", fill =NULL,
  title = "Perceived project benefits by partisanship") +
  scale_fill_party() +
  theme_classic(base_size = 10) +
  annotate("text", vjust =1, hjust=1.5, x = Inf, y = Inf, size =3, color = "grey40", label = paste0("N=",nrow(subset(g, !is.na(greenbenefit))))) +
  scale_y_continuous(expand = c(0, 0), limits = c(0, 1), labels = scales::percent_format(accuracy = 1)) +
  theme(
    legend.position = "inside",
    legend.position.inside = c(.15, .8)
    )
save(p.benefit.pid, file = here("data", "output", "fig_perceived_benefits_pid.rdata"))

# Diagnostic figure (not in paper)
if (isFALSE(pnas)) {
  save_pnas_pdf(p.benefit.pid, here("output", "pnas", "figures", "fig_perceived_benefits_pid.pdf"))
}

# Estimate models----

m.benefit <- list()
m.benefit[[1]] <- est_model(varname = "greenproj_bin+union_mem_z+eprice_z", outcome = "greenbenefit_bin", se = "state", fe = NULL, data.in = g)
m.benefit[[2]] <- est_model(varname = "greenproj_bin*pid3+union_mem_z+eprice_z", outcome = "greenbenefit_bin", se = "state", fe = NULL, data.in = g)
m.benefit[[3]] <- est_model(varname = "greenproj_bin*demshare_z+union_mem_z+eprice_z", outcome = "greenbenefit_bin", se = "state", fe = NULL, data.in = g)

benefit.cap <- paste0(
  "Linear probability model of perceived project benefits \\label{tab:perceived_benefits}"
)

benefit.notes <- paste0(
  "\\textit{Notes:} ",
  "Unit of analysis is the individual survey respondent. ",
  "Dependent variable = 1 if respondent perceives benefits from local green energy projects, 0 otherwise. ",
  "Estimates are OLS with cluster-robust standard errors by state. ",
  "Continuous covariates are standardized: county-level variables use within-state standardization; ",
  "state-level variables use full-sample z-scores; individual-level indices are standardized by construction. ",
  "* $p<0.05$, ** $p<0.01$, *** $p<0.001$."
)

tab.benefit <- here("output", "pnas", "tables", "tab_S27_perceived_benefits.tex")

modelsummary(
  m.benefit,
  title = benefit.cap,
  notes = benefit.notes,
  coef_map = c(
    coefmap,
    "greenproj_bin:pid3Neither" = "Visibility x Neither party",
    "greenproj_bin:pid3Republican" = "Visibility x Republican",
    "greenproj_bin:demshare_z" = "Visibility x 2020 county Biden share"
  ),
  stars = c("*" = 0.05, "**" = 0.01, "***" = 0.001),
  add_rows = data.frame(
    c("Sample Fixed Effects", "State Fixed Effects"),
    c("No", "No"),
    c("No", "No"),
    c("No", "No")
  ),
  gof_omit = "IC|RMSE|Std|FE|Within",
  gof_map = gm,
  output = "tinytable",
  label = "tab:perceived_benefits",
  escape = FALSE,
  fmt = fmt_significant(2)
) %>%
  theme_latex(resize_width = 0.38, resize_direction = "both", placement = "H") %>%
  save_tt(tab.benefit, overwrite = TRUE)

avg_slopes(m.benefit[[2]], variables = "greenproj_bin", by = "pid3")